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ABSTRACT 

A stochastic model for the energy of a fl aring solar active region i s pres ented, 



generalising and extending the approach of IWheatland Glukhovl (119981 ). The 
probability distribution for the free energy of an active region is described by 
the solution to a master equation involving deterministic energy input and ran- 
dom jump transitions downwards in energy (solar flares). It is shown how two 
observable distributions, the flare frequency-energy distribution and the flare 
waiting-time distribution, may be derived from the steady-state solution to the 
master equation, for given choices for the energy input and for the rates of flare 
transitions. An efficient method of numerical solution of the steady-state master 
equation is presented. Solutions appropriate for flaring, involving a constant rate 
of energy input and power-law distributed jump transition rates, are numerically 
investigated. The flare-like solutions exhibit power-law flare frequency-energy 
distributions below a high energy rollover, set by the largest energy the active 
region is likely to have. The solutions also exhibit approximately exponential 
(i.e. Poisson) waiting-time distributions, despite the rate of flaring depending on 
the free energy of the system. 

Subject headings: Sun: flares - - Sun: corona - - Sun: activity - - Methods: 
statistical 



Introduction 



Solar active regions are the sites of occurrence of most solar flares. Large active regions 
may persist for se yeral solar rotatio ns, and produce dozens of significant flares during a 
transit of the disk (jRichardsonl Il95ll ) . The largest flares are believed to involve the release 
of more than 10 27 J of stored magnetic energy, the energy appearing in accelerated particles, 
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heating, bulk motion of material, and radiation (jHudsonl Il99ll ). The process underlying 
flares is accepted to be magnetic reconnect ion, although the details of the process remain 
the subject of study (IPriest fc Forbes! |2002| ). 



The dynamical energy balance of solar active regions presents challenges to our under- 
standing. The energy liberated in flares is thought to be excess or 'free' magnetic energy 
associated with electric current systems in the solar atmospher e, but the origin of the cur 



rents, and hence the source of the energy, is not well understood (IMcClymont fc Fisher 



1987 



Tandberg-Hanssen & Emslidll988l : iLeka et al.lll996l : IWheatlandll2000al : ISchriiver et al.ll2005r ) . 
Two popular pictures for energy supply are first that existing coronal magnetic structures 
are twisted and sheared by photospheric motions, producing currents, and second that new 
current-carrying magnetic flux emerges through the photosphere. These mechanisms suggest 
that energy supply to an active region may be described as a continuous process, driven by 
slow photospheric and sub-photospheric motions. In contrast, flare energy release is rapid 
and unexpected. An important point is that the size of the downwards jump in energy 
associated with a flare may be very large, by comparison with the amount of stored active 
region energy. 

The understanding of active region energy balance is hampered by the inability to 
calculate coronal magnetic energy. In principle the magneto-hydrodynamic (MHD) virial 
theorem permits the calculation of magnetic energy from vector field values inferred in the 
chromosphere subject to the assumption that the field is everywhere force free (e.g. Metcalf 
et al. 2005). However few chromospheric vector field determinations are made, and the 
reliability of the method is unknown. Methods for modeling coronal magnetic fields from 



photospheric or chrom ospheric boundary conditions are being developed (IMetcalf et al.ll2008 



Schrijver et al.l 120081 ) . and these may permit estimation of coronal magnetic energy. There 
has also been recent progress in methods for estimating the rate of supply of energy to an 
active region directly from observations (e.g. Welsch et al. 2007). For example, in ideal MHD 
the Poynting flux is S = (— u x B) x B//io, where u is the fluid velocity and B is the magnetic 
field, and in principle this quantity may be estimated at the photosphere from observations. 

Although we lack detailed quantitative information about the rate of energy supply to, 
and the energy stored in, active regions, we do have detailed information about solar flare 
occurrence. The energy released in flares may be estimated (albeit subject to some error), 
and the rate of occurrence of flares is observed. Two related statistical properties of flares - 
the frequency-energy distribution, and the distribution of times between events - have been 
studied in some detail, as summarized below. 



Studies of th e frequency-energy distribution show that it is a 
decades in energy (Hudson 1991 ; Crosby. Aschwanden. &: Dennislll99~3 



jower law over many 



Aschwanden. Dennis. &: Benz 
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19981 ). Specifically the distribution may be written 



M{E) = AE-\ 



(1) 



where M{E) is the number of flares per unit time and per unit energy E, the factor A is a 
(time-dependent) measure of the total flaring rate, and 7 « 1.5. Typically this distribution 
is determined for flares from all active regions on t he Sun over some period of time, but it 
also appears to apply to individual active regions (IWheatlandl l2000bl ) . which suggests that 
the power law is intrinsic to the flare mechanism. A popular model explain ing the power 
law is the avalanche model (ILu fc Hamilton! 1 1 99 It ICharbonneau et al.ll200ll ). in which the 
magnetic field in the corona is assumed to be in a self-organized critical state, and subject 
to avalanches of small-scale reconnection events. 

There has been considerable interest in the flare waiting-time distribution, and more 
generally in waitin g-time distributions for what they re veal about underlying physics in a 
variety of systems (jSanchez. Newman fc Carrerad I2002T). Determinations of flare waiting 



time distributions have given varied results (IPearce. Rowe fc Yeune 



Wheatland. Sturrock fc McTiernanlll998l ; iBoffetta et al 



2001 



Moon et al. 



1999 



1993; 



Moon et al. 



Biesecker 



2001 



1994 



Wheatland 



200 



3). The results suggest that the observed distribution depends on the 



particular active re gion and on time , and that it is also influenced by event definition and se- 



lection procedures (IWheatland 



2001 



20051 ; iBaiesi. Paczuski. fc Stella 



2006 



Buchlin. Galtier. &: Vellill2005uPaczuski. Boettcher. &: Baiesi 
|). For some active regions, the distribution appears to 



be consiste nt with a simple P oisson process, i.e. independent events occurring at a constant 
mean rate (IMoon et al.l 120011 ). The corresponding waiting-time distribution is exponential. 
Other active regions show time-variation in the flaring process, and flare occurrence may 
be approximated by a piecewise constant, or more generally time- varying Poisson process 



( IWheatland 



200 1]) . On longer t ime scales a power-law tail is observed for events from the 



whole Sun faoffetta et al.lll999h. This may b e accounted for in terms of a time-dependent 



Poisson model (jwheatland &: Litvinenkq 2002), although so me authors have argued that the 



power law has fundamental significance (ILepreti et al.ll200ll ). We note that the waiting-time 
distribution has often been considered in isolation from other flare statistics, in particular 
the frequency-energy distribution. 

To understand the statistics of flare occurrence, it is desirable to have a general the- 
ory for the energetics of an active region, relating the free e nergy of the active regi on to 
the observed frequency-energy and waiting-time distributions. iRosner fc Vaianal (119781 ) pre- 
sented the first model of this kind, which was analogous to Fermi acceleration. In this model 
active regions experience exponential growth in free energy between flares which occur as 
a Poisson process in time, depleting all accumulated energy. This gives a power-law flare 
frequency-energy distribution for large energies, and presupposes Poisson occurrence. Litvi- 
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nenko (1994; 1996) generalized the model to incorporate different rates of supply of energy 
to the system, but the Litvinenko m odels retained the featur e that each flare releases all 
of the f ree energy. This aspect of the iRosner &: Vaianal (119781 ) model was criticized by iLu 
( 119951 ). IWheatland &: Glukhovl (119981 ) introduced a model permitting arbitrary changes in 
free energy at each flare. The model assumes the free energy E of an active region in- 
creases secularly between jump transitions downward in energy (flares), which occur at a 
rate a(E,E') for jumps from E to E', per unit energy. A master equation describes the 
steady-state energy balance, and the solution of this equation is the probability distribution 
P(E) for the free energy. The flare frequency-energy distribution is given by the convolution 



/"OO 

M{E) = / P(E')a(E', E' - E)dE'. 
Je 



(2) 



Wheatland &: Glukhovl (119981 ) investigated solutions to the master equation for constant en- 
ergy input and for a choice of transition rates a(E, E 1 ) ~ {E— -E") -7 , which leads to power-law 
behavior N(E) ~ E~^ (below a high energy rollo ver set by the highest energy an active region 
is likely to have). IWheatland &: Glukhovl (119981 ) did not determine a waiting-time distribu- 
tion for the model. It was argued that the derived power-law solutions are consistent with 
an avalanche-type model, and avalanche models have simple Poisson waiting-tim e statistics 

( Wheatland. Sturrock fc McTiernan 1998 ; Sanchez. Newman fc Carrerasl 2002 ^ However, 

this presents a puzzle: the total rate of flaring in the IWheatland fc Glukhovl (119981 ) model 
is given by 



X(E) 



a(E, E')dE', 



(3) 



which depends on the energy of the system. Hence it is expected that the occurrence of 
flares is not strictly Poisson, since the occurrence of a flare changes the energy of the system, 
and hence the instantaneous total rate of flaring. Non-Poisson waiting-time statistics might 
then be expected. 



Recently lDaly &: Porporatol (120071 ) demonstrated how to determine steady-state waiting- 



time d istributions for continuous time processes with arbitrary jump transitions. The lDaly fc Porporato 
(120071 ) theory is quite general, applying to any system described by a single time-dependent 
stochastic variable x{t) following a deterministic trajectory interrupted by positive or neg- 
ative jumps of random timing and size. The probability distribution P(x,t) for x(t) is de- 
scribed by a master equation, and the waiting-time distribution for jumps in the steady state 
may b e obtained from the solution to the steady-state master equation. iDaly fc Porporato 
(120071 ) demonstrated the theory in application to simple models for human attention, for volt- 
age across a nerve membrane, and for soil moisture content associated with rainfall events. 
In each case the models were analy t ic, iny olving simple solutions t o the master equation. 
In this paper the IDaly &: Porporato! (120071 ) theory is applied to the IWheatland fe Glukhov 
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( 119981 ) model for active region free energy, and waiting-time distributions are derived for so- 



lutions of relevance to flares. The application of the theory is relatively st raightforward: the 



Whea tland & Glukhov (1998) model is in the class of models considered by lDaly fc Porporato 



( 120071 ). the time-dependent stochastic variable being the active region energy E{t). Minor 
modifications to the theory are required because E is a positive definite quantity and because 
the jump transitions in E are always negative. A specific difficulty in applying the theory 
is that for the power-law form for rate transitions relevant for flares, the master equation is 
not amenable to analytic solution. In this paper an efficient numerical method of solution of 
the steady state master equation is presented and appli ed. This paper also cons i ders a more 



general form for flare transition rates than considered in IWheatland fc Glukhovl (119981 ). The 
results resolve the puzzle outlined above concerning whether the model produces an expo- 
nential waiting-time distribution. 

The sections of the paper are as follows. Section [2] presents the model, starting with the 
time-dependent master equation (§ 12. ip . Section [2T21 shows how steady-state waiting-time 
distribution may be obtained, § 12.31 presents simple analytic solutions illustrating the theory, 
and § 12.41 considers the information provided by moments of the master equation. Section [3] 
presents flare-like solutions to the master equation, starting with a justification of appropriate 
choices for the rates of transitions and for the energy supply rate (§ 13.11) . Section [3721 describes 
the numerical method, and § 13.31 presents the results. Section 0] discusses the results, and 
their significance for understanding solar flares. 



Model 



2.1. Master equation 



An active region is modeled as a system with free energy E(t) which evolves in time 
due to secular energy input and jumps downward in energy at random t imes and of random 



sizes. The syst em is described by the time-dependent master equation (IVan Kampen!ll992 



Gardinerl 120041 ) for the probability distribution P(E,t) of the free energy: 



dP(E,t) 
dt 



+ 



^[(3(E,t)P(E,t)]-X(E)P(E,t) 



P(E',t)a(E',E)dE', 



(4) 



where /3(E,t) describes the energy input rate at time t, a(E,E') describes the rate of flare 
jumps from E to E', and X(E) is the total rate of flaring, given by equation ([3]). The terms 
on the right hand side of equation (j3J) describe the system gradually increasing in energy 
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due to energy input, falling to a lower energy due to a flare, and falling from a higher energy 
d ue to a flare, respectively . This is the time-dependent version of the master equation given 



in 



Wheatland fc Glukhovl f|1998h . 



Following iDaly fc Porporatol (120071 ). we note that the system is also described by the 
stochastic differential equation 



dE 

— = P(E,t)-A(E,t) 



(5) 



where A(E,t) = Yli=i AEi8(t — U) describes accumulated losses in energy due to flaring, 
8(x) is the delta function, and where the times are given by a state-dependent Poisson 
process with occurrence rate X[E(t)]. The jump amplitudes AE are distributed according 
to the (state-dependent) distribution h(AE,E), defined by 



a(E, E - AE) = X(E)h{AE, E). 
so that J E h(AE, E)d(AE) = 1. 



(6) 



2.2. Steady-state waiting-time distribution 



Daly &: Porporatol (120071 ) showed how — assuming a steady state — the waiting-time 
distribution for the jump transitions may be derived. In this section we briefly re-iterate the 
theory, as it applies to the present model. 

Consider a deterministic trajectory described by equation ([5]), starting at energy E s and 
ending at a higher energy E e , the instant before a jump occurs. The distribution p e (E) of 
final energies E e is given by the rate of jumping at a given energy divided by the mean total 
rate of jumping, i.e. 

X(E)P(E) 

Pe{E) = 7^7 , (7) 

where P(E) is the steady-state solution to the master equation (TJJ and 

POO 

(A) = / X(E)P(E)dE (8) 
Jo 

is the mean total rate. The distribution p s (E) of starting energies E s is then given by p e {E) 
together with the distribution of jumps h(AE, E): 



POO 

p s (E)= / Pe (E')h(E'-E,E')dE' : 

J E 



(9) 
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which using equations (jSJ) and (j7]) may be rewritten as 

p s (£) = — / P(E')a(E',E)dE'. 
The waiting-time distribution is given by 



PAT) 



dT 

dr 



with 



where p\(E,t) is the solution to 



POD 

F{t)= / Px (E,r)dE : 
Jo 



9 -^§^- = ~ [P(E)p x (E,t)] - \{E)p x {E,t) 



(10) 

(11) 
(12) 

(13) 



with the initial condition p x (E, 0) = p s (E). Equation (fT3"|) describes the evolution of the 
system before a flaring jump occurs, i.e. over the deterministic trajectory starting at energy 
E s and ending at energy E e . 

A simpler form for the waiting-time distribution may be obtained when /3(E) = (3q, a 
constant. Solution of equation (TIB"]) by characteristics then gives 



Px(E,t) =p s (£- A,t)exp j-^ \ [E-Po(t-s)]d3 
assuming E > Pot, and p x (E,t) = otherwise. In this case 

Hr)- 

where 



f(u, t) = exp 
so the waiting-time distribution is 



p s (u)f(u,r)du 
\(/3qs + u)ds 



Pt{t) = / p s (u)\(I3qt + u)f(u, r)du. 



(14) 

(15) 
(16) 

(17) 



2.3. Steady-state analytic solutions 

Two analytic examples illustrate the application of the theory. The examples are not 
relevant for flares, because they do not produce power-law frequency-energy distributions, 
but they show how Poisson and non-Poisson waiting-time distributions may be obtained. 
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The case 0(E) = n a n d a(E , E') = cto (where ao and /3q are constants) was considered 



by IWheatland fc Glukhovl (119981 ). In this case equation (j3J) gives X(E) = a E, so the total 
rate of jumps is energy dependent and the waiting-time distribution will not correspond to 
a simple Poisson process. The analytic solution to the steady-state master equation is 

P{E) =aEe-^ aE \ (18) 

with a = a /p , and from equation (|2J) the frequency-energy distribution for jumps is a 
Gaussian: 

Af(E)=a e- 1 i aE2 . (19) 

\l/2 _____ 

From equation (jSj) the mean total rate is (A) = (a A)) > an d using equations (j7j) and ffTUl) 
the distributions of end- and start-energies for deterministic trajectories are 

p e (E) = f-Y^ a 3 / 2 E 2 e- l * aE2 (20) 



and 



Ps (E) =(-) aVV^ 2 (21) 



7T 



1/2 



respectively. Using equations (IT7|) and (J2T|) it follows that the waiting-time distribution is 
also a Gaussian: 

p T (r)= (^) 1/2 e-^oPor* (22) 



7T 



As a second example, we consider the case (3(E) = j3 and X(E) = A (where /3 and A 
are constants). Since the total rate is constant the waiting-time distribution must be 

p T (r) = A e- AoT , (23) 

i.e. jumps occur in time as a simple Poisson process. From equation (^ this case requires 
a(E,E') = Xo/E. The corresponding solution to the steady-state master equation is 

P(E) = b 2 Ee~ bE } (24) 

with b = Xo/Po, and from equation (|2J) the frequency-energy distribution of jumps is expo- 
nential: 

J\f(E) = bX e- bE . (25) 

The mean total rate of jumps given by equation (jSJ) is (A) = A , and using equations (j7j) 
and ffTOl) we have p e (E) = b 2 Ee~ bE and p s (E) = be~ bE . From equation (|17[) it follows that 
the waiting-time distribution is indeed given by equation (|23|) . 
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2.4. Moments of the master equation 

Mom ents of the master equation giv e useful information about the global behavior of 



solutions (IWheatland fc Litvinenkoll200ll ) 



The zeroth moment, obtained by integrating equation over all energies, gives the 
trivial result 

- j P(E,t)dE = 0, (26) 

i.e. normalization is preserved, provided (3(E,t)P(E,t) goes to zero as E — > and E — > oo. 

The first moment, obtained by multiplying equation by E and integrating over all 
energies, gives the simple statement of global energy balance: 

j t (E) = </?) - <r>, (27) 

where for any quantity / = f(E,t), the mean (/) is defined by 

POO 

(/>= / f(E,t)P(E,t)dE, (28) 

and where 



o 



r(E)= I (E — E')a(E, E')dE' (29) 
Jo 

is the total rate of energy loss at energy E. Equation (127|) requires (3(E,t)P(E,t) to go to 
zero as E — > and E — > oo. In the steady state equation (j27|) gives 

(/?) = (r). (30) 



3. Flare-like solutions 



In the following we consider solutions to the steady-state master equation which may 
be of relevance to solar flares. 



3.1. Choices appropriate for flares 

We restrict attention to the case (3(E) — /3q, a constant. The motivation is that active 
regions are externally driven, i.e. the energy is supplied from the sub-photosphere by external 
processes. In the absence of a back reaction, it is then expected that the energy supply rate 
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does not depend on the state of the system. In passing we note that in general the energy 
supply rate may depend on time. However, in this section we consider only steady-state 
solutions to the master equation. We return to this point in § 4. 



We consider the form 



a(E, E') = a E 5 (E - E'Y 1 6{E - E' - E c ) 



(31) 



for the flare transition rate, where E c is a low-energy cutoff , and 8(x) is the step function. The 
case 5 = was considered in I Wheatland fc Glukhovl (119981 ). The motivation for equation (I3T1) 
is that it may describe an avalanche-type system, in which energy transitions are intrinsically 
power-law distributed. The power law is assumed to originate in the microphysics of the flare 
process, and must be assumed at the level of this model. [Thi s is in contrast to models which 
attempt to account for the power law, e.g. iRosner fc Vaianal (119781 ).] The low-energy cutoff 
E c is needed to ensure X(E) is finite. The E s factor represents a possible dependence of the 
transition rate on the energy of the system. It is plausible that an avalanche-type system 
with more energy is more likely to contain unstable sites, and hence will flare at a higher 
rate. In the following we take 7 = 1. 5 in every instance, and consider two cases: 5 = 0, 



following IWheatland fc Glukhovl (119981 ); and 5 = 1. The choice of transition rates (13T|) leads 
to the total flaring rate 

X(E) = ^jE 5 (E^' +1 - E-i +1 ) . (32) 

Hence the rate of occurrence of flares is energy-dependent, and the waiting-time distribution 
will not correspond to a simple Poisson process, as pointed out in § 1. 



20021). Specifically, 



Application of the first moment with the choice of transition rates (13T1) lea ds to a simple 
estimate for the mean energy of the system (IWheatland fc Litvinenkol 
from equation (T29l) we have 

"0 _^5+2- 7 



r(E) 



2-7 



(33) 



for E 3> E c . Taking averages and making the approximation (E s+2 7 ) « (E) s+2 7 together 
with (P) = (3q in equation (1301) gives 



(E) 



7 



<5+2-7 



(34) 



Substituting equation (13~T1) into equation (j2J) leads to the form for the flare frequency- 
energy distribution 

POO 

M(E) = a Q E~~ 1 (E'f P(E')dE' (35) 

J E 
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for E > E c . Hence it follows that the frequency-energy distribution will be a power law 
with index 7 up to energies E at which P(E) becomes very small. Equation (1341) provides a 
crude estimate (a lower bound) for the energy at which the frequency-energy distribution is 
expected to depart from power law behavior. 



3.2. Numerical method 

The steady-state master equation may be non-dimensionalized by introducing new vari- 
ables E = E/Eq, P = PE , t = Pot/Eo, and a = E^a/fto, where E is a chosen scale for 
energy. [For the solutions corresponding to equation (13"T1) we take E = E c ] This procedure 
gives 

dP 

Je 

where 

A= / a(E,E')dE'. (37) 
Jo 

Hereafter we assume non-dimensional equations, but omit the bars. 



POO 

+ XP- _ P(E')a(E',E)dE' = 0, (36) 
Je 



E 



Equation (1361) is linear in P(E), and hence may be solved by discretizing in energy and 
solving a coupled system of linear equations. These equations must be supplemented by the 
normalization condition on P(E). Direct back-substitution provide s an efficient method of 



solutio n, and the details of the procedure are given in Appendix A. In I Wheatland fc Glukhov 



(119981 ) the steady-state master equation was solved by a relaxation procedure, but the ap- 



proach given here is more numerically efficient. 

The flare frequency-energy distribution is obtained from the solution for P(E) via nu- 
merical evaluation of equation §2§. The waiting-time distribution is similarly determined 
via numerical evaluation of equation (I17p . using an analytic form for f(u,r) obtained from 
equation (|T6|) . All numerical integrations use the extended trapezoidal rule. The numerical 
solution was tested on the analytic cases given in § 12.31 



3.3. Results 



First we consider the case 5 = 0, following IWheatland fc Glukhovl (119981 ) . Figured] 
illustrates the numerical solution of the steady-s tate master equation (|36l) for t he case «o = 
0.1, which is one of the two cases considered in IWheatland fc Glukhovl (119981 ). The upper 
panel shows the probability distribution P(E) for active region energy (as a linear-log plot), 
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the middle panel shows the flare frequency-energy distribution J\f(E) (as a log- log plot), 
and the lower panel shows the waiting-time distribution p T (r) (as a log-linear plot). As 
explained in § 13.11 the frequency-energy distribution is expected to be a power law with index 
7 = 1.5 below energies at which P(E) becomes small, and the expression (J3ID provides a 
lower bound for the departure from power-law behavior. The lower bound is shown in 
the uppe r and middle panels by a ver tical line. The upper and middle panels confirm the 



results of lWheatland fc Glukhovl (119981 ). The lower panel shows the waiting-time distribution 
Pt(t) (solid curve) as well as the Poisson distribution (A)e _ ^ r (dotted line) corresponding 
to the mean rate of flaring implied by the form of \(E) and the solution for P(E) [see 
equation (|SJ)]. Note that the units for time in the lower panel are E c / /3q, following the non- 
dimensionalization in § 13.21 The waiting-time distribution for the model is approximately 
Poisson, although there is a slight deficiency of long waiting-times by comparison with the 
Poisson distribution. 

Figure [2] shows the c a se S = and «o = 0.02, which is the other case considered 



by IWheatland fc Glukhovl (119981 ). in the same format as Figure [H For a lower value of 
a the system flares less often and hence is more likely to have larger energy. Hence the 
distribution P(E) shown in the upper panel is shifted to higher energy and has a higher 
mean. The flare frequency-energy distribution (middle panel) is a power law over more 
de cades in energy than fo r the c ase ao = 0.1. These results are consistent with the findings 



of lWheatland fc Glukhovl (119981 ) . The lower panel shows the waiting-time distribution (solid 
curve) as well as the Poisson distribution corresponding to (A) (dotted), although the two 
curves are almost indistinguishable. 

The results in Figures [1] and [2] suggest that the 5 = model has a waiting-time dis- 
tribution which is close to being strictly Poisson (exponential), and that the approximation 
becomes better for smaller values of otQ. This may be understood in terms of the expres- 
sion (1321) . For E ^> E c the total rate is X(E) ~ olqE~ 1+1 / (7 — 1), which is constant, in which 
case a Poisson waiting-time distribution is expected. For smaller values of ao the system 
is more likely to have larger energy, and hence the approximation E ^> E c will be better. 
Figure [3] illustrates this explanation for the case 6 = 0, a = 0.02. The solid curve shows the 
total rate as a function of energy, and the dashed line shows the mean total rate (A). The 
energy distribution P(E) is shown, with arbitrary normalization, by the dotted curve. We 
see that, over the range of energy for which the distribution P(E) is substantial, the total 
rate is approximately constant and equal to the mean total rate. These results resolve the 
puzzle identified in §[TJ the waiting-time distribution for the model is not strictly Poisson, 
but is a good approximation to an exponential. 

Next we consider the case 5 = 1, to examine what happens when the rate of flare 
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transitions increases with the energy of the system. Figures H] and [5] show the cases ao = 10~ 3 
and otQ = 10~ 5 respectively, in the same format as Figures 1 and 2. First consider Figure HI 
The energy distribution P(E) shown in the upper panel is qualitatively similar to the 5 = 
case, although the distribution declines more rapidly at higher energies, so that it is more 
skewed in a linear-log representation. Since the rate of transitions increases with energy, 
the system is less likely to be found at very large energies, and this explains the rapid 
decline. The middle panel shows the flare frequency-energy distribution Af(E), which is a 
power law with index 7 below a high energy rollover. The estimate f)34p for the mean of the 
distribution (vertical line) provides a lower bound for the departure from power law behavior. 
The lower panel shows the waiting-time distribution p t (t), and the Poisson distribution 
implied by (A). The waiting-time distribution is approximately Poisson, but has a deficit 
of long waiting times. Figure illustrates the case with reduced flare transition rates. The 
distribution P(E) (upper panel) is shifted to higher energies, and is again quite skewed in 
the linear-log representation. The frequency-energy distribution Af(E) (middle panel) is a 
power law over more decades in energy. The waiting-time distribution p T (r) (lower panel) is 
again approximately exponential, but departs somewhat from the Poisson model, including 
showing an excess of long waiting times. 

The results in the lower panels of Figures @] and [5] suggest that, for the 5 = 1 model, the 
waiting-time distribution is approximately exponential but shows some departure from the 
Poisson case depending on the parameters of the model. The approximate Poisson behavior 
is perhaps surprising because in this case the total rate of flaring [given by equation (|32p ] 
varies approximately linearly with E (for E ^> E c ) . Figure M illustrates this for the case 
5 = 1, ao = 1CT 5 , using the same format as Figure [31 The total rate (solid curve) increases 
substantially over the range of energies the system is likely to have [the dotted curve shows 
P(E)}, and may be substantially different to (A) (the dashed line). However, p T {j) is defined 
by a complicated average of the rate over energy, which is different for different waiting 
times [see equation (IT7j) ]. and the numerical results show that the resulting waiting-time 
distribution is approximately Poisson. 

4. Discussion 

A stochastic model is presented for the free magnetic energy of a flaring solar active 
region. The energy of an active region is assumed to grow deterministically between random 
flare events at which the energy jumps downwards by an amount equal to the flare energy. 
Flares jumps occur from energy E to E' with a rate a(E, E') per unit time and per unit 
energy, and energy input occurs at a rate (3{E). Active region energy is then described by 
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a distribution P(E) which is the steady state solution to a master equation. This distri- 
bution determines two observable distributions, namely the flare frequency-energy distribu- 
tion, and the waiting-time dis tribution. The model generalizes and extends the approach of 



Wheatland fc Glukhovl (119981 ). Novel aspects of the work presented here i nclude the determi- 



nation of waiting-time distributions [following general theory presented by lDaly fc Porporato 
(120071 )]. consideration of a more general form for the rate of flare transitions, and introduction 
of an efficient method of numerical solution of the steady-state master equation. 

The form a(E,E') = o,qE 5 (E — E')" 1 ' 5 for flare tra nsitions is investigated, for th e 
cases 5 = and 5 = 1. The case 5 = was considered by IWheatland & Glukhovl (119981 ). 
motivated by the avalanche model. For both cases the model is shown to produce power- 
law flare frequency-energy distributions below a rollover at high energies, due to the active 
region having a finite energy. For both cases the waiting-time distribution is approximately 
exponential (Poisson). For the case 5 = this may be understood in that the total rate 
\{E) is approximately constant for E 3> E c , which becomes a good approximation for small 
«o, when the system is more likely to be found at large energies. This result is consistent 
with the interpretation of this model as avalanche-like, since avalanche models have simple 
Poisson statistics. For the case 5 = 1 the interpretation is more complicated, because the 
total rate varies approximately linearly with energy. However, the waiting-time distribution 
is determined by an average of rates over the possible energies of the active region, and 
numerical evaluation shows that the result is approximately exponential. 

The general model introduced in § 12. II includes time dependence in the driving rate, but 
we have focused on the steady state throughout this paper. Ma ny active regions exhibit 
large variations in the rate of flaring during a transit of the disk (IWheatlandl 120011) . This 
beha vior is often linked e.g. with the emergence of new magnetic flux (IRomano fc Zuccarello 



20071 ). which suggest that it is a response to an increased rate of driving. Hence we have 



neglected an important aspect of active region energetics. Time-dependent driving will 
influence the observed waiting-time distribution. In the simplest case time variation might 
be represented by a piecewise constant variation in the driving rate. If the system adjusts 
suitably quickly to changes in driving, the steady state solution applies to each piece. The 
waiting-time distribution is then a weighted s um over the steady-state distributions applying 
to each piece (IWheatland fc Litvinenkoll200ll ). Based on the results presented in this paper, 
the waiting-time distribution for active regions is expected to be approximately exponential 
provided the rate of driving of the system is constant. If the rate of driving is time varying, 
the distribution will depart from exponential. A time-dependent model will be investigated 
in future work. 



One shortcoming of the model is that it does not describe energy loss from the system 
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by mechanisms other than flaring, for example loss due to flux submergence, or quasi-steady 
background dissipation. However, a simple generalization of the master equation permits 
this. Specifically, if the secular energy 'input' is replaced by small jumps in energy (which 
may be positive or negative), then the energy gains and losses may be represented by Fokker- 
Planck terms in a generalized master or Chapman-Kolmogorov equation. Specifically, the 
secular energy increase term —d [j3(E,t)P(E,t)] /dE on the right hand side of equation (jlj) 
may be replaced by a pair of terms -d [(3 1 (E,t)P(E,t)} /dE + \d 2 [f3 2 (E,t)P(E,t)\ /dE 2 , 
where the coefficients f3i(E,t), i = 1, 2 repr esent first and second moments of energy changes 
associated with the small jump transitions (IVan Kampenj|l992l ; lGardinerll2004j ). It is straight- 
forward to solve the resulting generalize master equation by discretization and solution of 
the resultin g linear system, extending the approach presented in the Appendix. However, in 
this case the lDaly fc Porporatol (120071 ) method for determining the steady-state waiting-time 
distribution needs to be modified. This model will be investigated more completely in future 
work. 

The results presented here show how it is possible to construct a model for active region 
energetics which directly predicts observable flare statistics, namely the flare frequency- 
energy, and waiting-time distributions. In principle the observations may be used to deter- 
mine the energy supply and flare energy release terms in the model, i.e. f3(E, t) and a(E, E'). 
However, the observations are not really precise or unambiguous enough to identify these 
terms with certainty. In particular, the interpretation of the waiting-time distribution is 
complicated by the time dependence of the energy supply. The situation would be improved 
by an ability to estimate coronal free energy and the rate of supply of energy to active re- 
gions from observ ations. Reliable methods for estimating these quantities are the subject of 
current research (IWelsch et al.l 120071 ; ISchrijver et al.ll2008l ). If such methods are developed, 
the theory developed here will be of greater significance. It may provide valuable insight 
into the flare mechanism, as well as being of practical benefit for flare prediction. 



The author thanks Ian Craig for pointing out that the master equation can be solved 
as a linear system by back-substitution, and a referee for helpful comments which improved 
the presentation. 
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A. Appendix 
Numerical solution of the steady-state master equation 

Discretizing equation (155]) at energies E t = iA gives 

^ N-2 

+ XiPi -~ A Y1 ( p i a it + p J+i a i+i.O = ( A1 ) 

where 



p p .N-2 

2A ' 2' 



A. = ~ A + (A2) 

3=0 

and where P, = P(P«) and a^- = a(Ei,Ej). Centered differencing is used for the deriva- 
tive, and the extended trapezoidal rule is used for the integrals. Equation (1A1I) with 
i = 1, 2, N — 1 (and the assumption P/v = 0) may be supplemented by the normalization 
condition 

N-2 

-A ( p + p +i) = 1 (A3) 



2 

i=0 



to give N linear equations in the N unknowns Pq, Pi, P2, P, 
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The resulting linear system may be solved efficiently by back-substitution as follows. 
Equation ( 1A1I) may be re-written 



Pi-i = P l+ i + 2A 



^ AT-2 

AjPj - -A (PjQtj.i + Pj+1%+1, 



2 

3=* 



(A4) 



which expresses Pj_i in terms of Pj, Pj+i, Pj+2, •••> P/v-2- Hence if we assume a value for P/v-i, 
we can apply equation (IA4I) to solve for P/v-2, and then apply it again to solve for P/v-3, 
etc. In this way we can determine P/v-i, P/v-2, P/v-3, Po up to an unknown normalization 
factor. The factor is determined by applying equation (1A3j) . Specifically, the solution is 
given by P/ (with z = 0, 1, 2, ...,N — 1), where 

Pi 

* = Jae^^ + Ph-i)" (A5) 
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Fig. 1. — Numerical solution to the steady-state m aster equation for the case S — 0, 7 = 
1.5, and «o = 0.1, one of the cases considered in IWheatland fc Glukhovl (119981 ). Upper 
panel: probability distribution for free energy P{E); middle panel: flare frequency-energy 
distribution J\f(E); lower panel: flare waiting-time distribution p T (r). 
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Fig. 2. — Numerical solution to the steady-state m aster equation for the case S = 0, 7 = 
1.5, and «o = 0.02, one of the cases considered in IWheatland fc Glukhovl (119981 ). Upper 
panel: probability distribution for free energy P{E); middle panel: flare frequency-energy 
distribution J\f(E); lower panel: flare waiting-time distribution p T (r). 
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Fig. 3. — The total rate of flaring X(E) versus energy (solid curve) for the case 5 = 0, 
7 = 1.5, and a = 0.02, and the mean total rate (A) (dashed line). The energy distribution 
P(E) is also shown, with an arbitrary normalization (dotted curve). 
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Fig. 4. — Numerical solution to the steady-state master equation for the case 5 — 1, 7 = 1.5, 
and ao = 10~ 3 . Upper panel: probability distribution for free energy P(E); middle panel: 
flare frequency-energy distribution Af(E); lower panel: flare waiting-time distribution p T (r). 



-23- 




10000 



Active region energy 



10° 



Ld 




1 10 100 1000 10000 

Flare energy 



qj 10" 1 




o 



i m- 6 l . . . .. I 

100 200 300 400 

Waiting time 



Fig. 5. — Numerical solution to the steady-state master equation for the case 5 — 1, 7 = 1.5, 
and ao = 10~ 5 . Upper panel: probability distribution for free energy P(E); middle panel: 
flare frequency-energy distribution Af(E); lower panel: flare waiting-time distribution p T (r). 
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Fig. 6. — The total rate of flaring X(E) versus energy (solid curve) for the case 5 — 1, 
7 = 1.5, and a = 0.02, and the mean total rate (A) (dashed line). The energy distribution 
P(E) is also shown, with an arbitrary normalization (dotted curve). 



